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ABSTRACT 



This thesis proposes a method of identifying the dynamics of ship 
angular motion at sea as the basis for ship motion estimation. Various 
sources of information are discussed. A digital model to simulate 
ship's motion is developed. Identification algorithms are presented 
with the results of their stochastic digital simulation. Sensitivity of 
identification error to sample rate was investigated. A plan for ship- 
board implementation utilizing a digital computer is presented. 
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INTRODUCTION 

The sophistication of modern shipboard systems is imposing 
accuracy requirements in measuring a ship's attitude. A 10 mil error 
in yaw correction to a gun platform is a 200 yard error at 10 nautical 
miles. System time delays and additive measurement noise are key 
sources of system error. Automatic carrier landing systems require 
that the predicted touch down point on the deck be known. Present 
system effectiveness degenerates rapidly as sea states increase. 

When the wave-off decision point is reached by the aircraft, the auto- 
mated landing system must know the aircraft position relative to the 
desired impact point at the time of touchdown. If pitch and yaw can 
accurately be predicted, a true all weather capability may be antici- 
pated. For these and many other applications the accurate knowledge 
of present position and the prediction of future positions are key 
elements of the system environment. 

An accurate prediction of ship motion must necessarily involve 
information from the forcing function. Proper interpretation should lead 
to an identification of the current dynamic structure of the ocean waves. 
This could yield a more accurate input to the prediction of sea states 
and ocean weather propagation. 

The French have done considerable work in the field of spectro- 
angular ocean wave analysis [1,2]. However their work is success- 
ful only in predicting mean and extreme wave amplitudes and periods. 
The primary available inputs are local surface wind conditions dis- 
tributed over the ocean area. Work in this area is not yet sufficiently 
advanced to be applied to a real time estimation problem. 

The investigation of the statistical properties of ocean waves has 
been a well traveled route. R. L. Wiegel [ 3 ] from the University 
of California has done extensive work using both surface and bottom 
moored pressure sensing units. Data reduction again shows normal 
distributions of wave period and amplitude in open waters. 
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Ship motion estimation may be approached as a stochastic con- 
trol problem. Thus the ship is the plant and the ocean wave is the 
forcing function with normal Gaussian distribution. Pitch, roll, 
and yaw will be treated as separate, uncoupled subsystems dis- 
cretely sampled to facilitate digital simulation and data processing. 
Only the roll mode will be referred to in the remainder of this paper 
on the assumption that the other modes may be estimated by the same 
techniques . 
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PLANT SIMULATION 

To evaluate a technique for ship motion estimation, it is first 

necessary to have a source of data. A fourth order plant composed of 

two cascaded identical second order systems was selected for motion 

simulation. This plant has known natural frequency ( a? ) and damping 

n 

coefficient ( C ), 
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Vs + 2 C W n S + 00 n J 
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determines the dominant characteristics. 



n 

discrete recursive state equation in matrix form is 



X(k + 1) = 0 X(k) + T W(k) . 



( 1 ) 



( 2 ) 



The 

(3) 



W(k) is the discrete, white Gaussian excitation with zero mean. The 
desired system output is the discrete time sequence of roll angle 
corresponding to X^ where 

X 1 (k) = HX(k) , H = (1 0 0 ... ) . 

Additive measurement noise, v(k), may be included at the output to 
complete the simulation model, fig. 1. V(k) is assumed to be white 
Gaussian noise, uncorrelated with W(k). The process is described by 
z(k) = HX(k) + v(k). (4) 

It will be convenient to transform 0 into the canonic companion 
matrix form, 0 *, using Lee's observability criteria [4 ] , giving 
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H 



H 0 



D = 



H 0 



n-1 
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0 ! 



I 



0 * — D 0 D 



L 



(5) 




and the companion distribution matrix 



r* = Dr. 



(6) 



From the work of Lee and Ho [ 5 ] this defines a new recursion equation 
with the same output values as equation (4) as seen below 



where W(k), H, v(k), and z(k) are the same but the Y elements are 
different. 

Representing the simulation model in this form reduces both the 
number of computations and the time required to generate a sample set 
of data on a digital computer. Since the data set is unchanged by the 
foregoing manipulations, the dominant characteristics are preserved. 
Values of f= 0.7 and u.' =0.2 tt were used as reasonable estimates 



for the roll dynamics of a moderate size ship. The excitation, W(k), 
was a random number generating subroutine with zero mean and stan- 
dard deviation of one. Multiplying the random numbers by a con- 
stant permitted simulation of any desired sea state. An example of 
roll angles generated using a constant of twenty and a sample period 
of 0.1 second is shown in fig. 2. 



Y (k+ 1 ) = 0 * Y(k) + T * W(k) 



z(k) = H Y(k) + v(k) 



(7) 



n 
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Fig. 1. Block Diagram of the Simulation Model. 




Fig. 2. Roll Angles vs Time (10 sec/in) Generated by the Simulation 
Model Excited with Discrete Gaussian Noise, a= 20, at the 
Sampling Interval, 0.1 sec. Dominant Period = 10 sec, 
r = 0.7. 
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IDENTIFICATION 

The problem of identification arises when the set of differential 
equations describing the system dynamics are not known, uncertain, 
or are subject to periodic or continuous change. On a large time 
scale the ship dynamics may be expected to change as fuel, stores, 
ammunition, etc. vary in quantity and location. Short time changes 
may be anticipated due to changes in ships heading and speed. 

In state vector notation the recursion equation for ship roll 
motion (uncoupled) will have the form; 

X(k+1) = $X(k) + rw(k) (8) 

where W(k) is the discrete Gaussian excitation. $ is the fixed 
state transition matrix of rank N, the order of the system. The term 
fixed is used here on the assumption that the changes in dynamics 
mentioned above occur at such a sufficiently slow rate that $ may be 
adequately represented as time invariant over a finite time interval. 
Thus $ must be periodically identified to follow ship motion within 
accuracy specifications. 
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IDENTIFICATION OF THE FREE DYNAMIC SYSTEM 
A free dynamic system may be represented by the equation 

Y(k+1) = 0 * Y(k) . (9) 



Lee developed an identification scheme for a plant having no numerator 
dynamics [4 ] . Lee's development gave the following equations; 
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and z(k) = HY(k), a scalar. 

Using the same data, the order of the system may also be iden- 
tified. Choose some M greater than N and build an M x L matrix, A. 
L = 1,2,3,.. .N,N+1, 



A = 



Z 1 Z 2 Z £ 

Z 2 Z 3 Z £ +1 

* • * 

• • * 

Z m Z m+1 Z m+l - 1 



T 

The product A A will be positive definite for L less than or equal to N 
and singular for L greater than N. Hence the order of the system is 
the maximum nonsingular value of L. 

A problem arises from the possible occurrence of numerator dy- 
namics which leads to residues (errors) in the identification. In z- 

th 

transform theory the N row of 0 * corresponds to the denominator 
coefficients of the system z-transform. The numerator remains 



13 



unidentified. Hence, if there is an unknown non- scalar numerator in 
the real plant corresponding to a nonzero initial condition, errors will 
result. 

Identification using this method was investigated with excellent 
results by setting an initial condition for roll angle and then releasing 
the plant. If excitation is applied while measurements are being taken, 
the plant is no longer a free dynamic system. Identification attempted 
with random excitation applied was inaccurate and unreliable. 
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IDENTIFICATION FOR THE STOCHASTIC CASE 
Using the transformed recursion equation in the new state space, 
Lt. Ralph Hudson's unpublished doctoral notes show the development 
of a straight forward identification scheme as follows: 



Y(k+1) = 0 * Y(k) + co(k) 



(ID 



where 



Y(k) = 



"k-n+1 



'k-1 



, co (k) = T *W(k) , 



and the scalar, W(k), represents the white Gaussian excitation. If 

T 

equation (11) is post multiplied by Y(k) ; 



Y(k+l)Y(k) T = 0 *Y(k)Y(k) T + co (k)Y(k) T 



( 12 ) 



Taking the expectation of both sides; 

E[Y(k+l)Y(k) T ] = 0 *E[Y(k)Y(k) T ] +E[co(k)Y(k) T ] (13) 

T 

Notice that the E[ Y(k+l)Y(k) ] is just the autocorrelation function, 

T 

R( t), for tau equal to one. Similarly, E[Y(k)Y(k) ] is the auto- 
correlation for tau equal to zero. Therefore 0 * = R(1)R(0) ^ . 

If the order of the identification is M and M £N, the true system 

th 

order, the inverse will exist. This is implied by the fact that an N 
order linear system may be defined by N linear independent differ- 
ential equations. Hence any M ^ N of these equations are also in- 
dependent. If M >N independence is lost and R(0) of rank M >N is 
singular. Thus if N is unknown, it may be identified by testing for 
the largest non-singular rank of R(0). 

For application to discrete linear (Kalman) filtering it is also 
useful to note that the Q required for determining the gain matrix may 
be concurrently identified as follows; 

Q = T * E[ W(k)W(k) T ] r * T = E [ co (k) co (k) T ] (14) 
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Hence 



Q = E [ { Y(k+1) - 0 *Y(k) } { Y(k+1) - 0*Y(k)} T ] 
which simplifies to 

Q = R(0) - 0 *R(1) T . 



(15) 



This scheme was programmed in Fortran 63 for simulation on the 
CDC 1604 digital computer using data from the simulation model. 

Batch processing and recursive identification were investigated. 

1. Batch processing: A sample set of roll angles was generated 

and stored. The stored data was then batch processed to form the auto- 
correlation functions, R(0) and R(l), for the set. A matrix inversion 
routine was used to invert R(0). 0* was then identified as the matrix 

product, R(1)R(0) 

If r * was restricted to a single non zero element, identifi- 
cation of the fourth order plant was good using 50 samples and ex- 
cellent for 500 samples. However, if T * had all non zero elements, 
residues due to numerator dynamics caused errors as great as 110% 
of the true values of the elements. Increasing the sample size from 
500 to 900 showed no improvement. 

2. Recursive identification: Hudson's identification scheme may 
be manipulated into the following recursion equations by using a matrix 
inversion lemma. 



0 *(k+l) = 0*(k) + (y(k) T P(k)Y(k) + ” 1 ^Y(k+1) - 0 A *(k)Y(k)^ Y(k) T P(k) 



P(k+1) = P(k) 



^1 - [ Y(k) T P(k)Y(k) + 1] 1 Y(k)Y(k) T P(k) 



where P(k) = R(0) 



-1 



Let 0 (k) = Y(k) P(k)Y(k) + 1 



)"■ 



a scalar. 



Then 

0 A *(k+l) = 0 A *(k) + 0 (k) (yOc+ 1) - 0*(k)Y(k)) Y(k) T P(k) 



(16) 
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and 



P(k+1) = P(k) 



0 - 



jS (k)Y(k)Y(k) P(k) 



) 



(17) 



2 

The above relationships identify an N by N 0 * matrix of N elements. 

For the scalar observable case N(N-l) of these elements are already 

2 

known to be zeros or ones . Hence N elements are identified to learn 

th 

N unknowns in the N row. Thus a more efficient method of the fore- 
going has been derived by Yu-Chi Ho [ 5 ] and R.C.K. Lee [ 4,5 ] 

th. 

which recursively estimates the elements of the N row. 

0 (k+1) = 0 (k) + P(k)Y(k) { z (k+1 ) - 0(k) T Y(k)j 0 (k) (18) 



P(k) = P(k+1) 



( 



I - 0 



(k-l)Y(k-l)Y(k-l) P(k 



■«) 



(19) 



where z(k+l) = the scalar measurement at stage (k+1) 

a a T th 

0 (k) = (Nxl) column such that 0 (k) is equal to the N 

row of 0 *(k) 

Initialization of the recursive equations might be accomplished for 
a -i 

0*(O) = S S or by using a small sample set of 2N+1 or more 

measurements and forming R(l) and R(0) where by 0*(O) = R(1)R(0) 
and P(0) = R(0) ^ . Both of these methods require the complicated matrix 
inversion routine. From Lee's work [4 ] the convergence rate appears 

A A 

to be satisfactory for any reasonable estimate of 0*(O) or 0 (0) . Ex- 
perience seems to bear this out for the simulation model. The P matrix 

A 

ought to be inversely proportional to time approaching zero as 0 
approaches 0 . P has its most significant effect on the initial rate of 
convergence. Experience indicates that if P is initialized with a 

g 

reasonably large set of values on the main diagonal (~10 ), con- 
vergence is satisfactory. Additional work in this area would be useful 
to define an optimal P(0), perhaps in terms of initial convergence rate, 
measurement noise statistics, and identification error as a function of 
time. 

Convergence tests were run using the scalar observable roll angles 
generated by the fourth order model. A ten second period, damping 
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coefficient of 0.7 and an excitation constant of twenty were used. The 
sample period was one second. Three cases were investigated using 
the same sample set. 

The first case was with T * suppressed to zeros in all but the last 
element to eliminate numerator dynamics . The convergence of 600 
samples is ’shown in fig. 3a,b,c,d for each of the fourth row elements 
of 0 *. The second case was the same as the first but with Gaussian 
additive measurement noise having a standard deviation of one. Ex- 
cellent convergence is shown in fig. 4a,b,c,d. The third case used 
the entire T * vector and no measurement noise. The residues are 
apparent in fig. 5a,b,c,d. Comparing the eigenvalues of 0 * with 

A 

those of 0* with residues indicates a considerable shift, fig. 6. The 
responses to different sample rates in the following section indicates 
that the residues may be minimized sufficiently to provide acceptable 
root locations . 

The true computed values for the plant with a 10 sec period sam- 
pled at 1 sec intervals are as follows. 
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Fig. 3. Recursive Convergence of $* vs Sample Number, No Measurement 
Noise, 10 sec Dominant Period Plant with Zero Initial Conditions. 
Discrete Gaussian Excitation, < 7 = 20 , Applied at the Sample 
Interval of 1 sec . F * Suppressed (r*l,2,3 = 0),F = 0.7. 

The Straight Line Represents the True Value of the Particular 
0 * Element, (a) £ * (4 , 1) vs Sample Number . (b) $*(4,2) vs 
Sample Number, (c) 0 *(4,3) vs Sample Number, (d) 0 *(4,4) 
vs Sample Number. 

(a) £*(4,1) vs Sample Number, (0 *(4, 1) =-0.17217) 
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(d) 0*(4,4) vs Sample Number, (0 *(4,4) = 2.3215) 
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Fig. 4. Recursive Convergence of ti * vs Sample Number with Additive 

Gaussian Measurement Noise, a = 1, using the 10 sec Dominant 
Period Plant, Zero Initial Conditions. Sampled and Excited at 
1 sec Intervals , T * Suppressed ( r*l , 2 , 3 = 0) , r = 0 . 7 . The 
Straight Line Represents the True Value of the Particular 0 * 
Element, (a) 0 *(4,1) vs Sample Number, (b) 0 *(4,2) vs Sample 
Number, (c) 0 *(4,3) vs Sample Number, (d) 0 * (4 , 4) vs Sample 
Number . 

(a) % *(4,1) vs Sample Number, (0*(4,1) = -0.17217) 
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(b) % *(4,2) vs Sample Number, (0 *(4,2) = 0.96328) 



24 






001 



002 



003 



004 



005 



X-SCflLE - J.00E+0Q UNITS/INCH. 
Y-aCflLE - 5.00E-01 UNITS /INCH. 



I 




25 



sic cio 



§ 



S -SCALE - 1.00E+82 UWT3/MH 
Y -SCALE - 5.00E-01 UMITS/IHCH 



t 







i 

§ 

000 801 00Q 003 004 006 



(d) 0 *(4,4) vs Sample Number, (0 *(4,4) = 2.3215) 
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Fig. 5. Recursive Convergence of 0 * vs Sample Number, No Noise 
using the 10 sec Dominant Period Plant, Zero Initial Con- 
ditions. Sampled and Excited at 1 sec Intervals, l" = 0.7. 

The Straight Line Represents the True Value of the Particular 
0 * Element. The Full T * vector was used to Demonstrate 
the Residues from Numerator Dynamics, (a) 0 *(4,1) vs Sample 
Number, (b) 0 *(4,2) vs Sample Number, (c) $*(4,3) vs Sample 
Number, (d) $*(4,4) vs Sample Number. 

(a) 0 *(4,1) vs Sample Number, (0 *(4,1) = -0.17217) 
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(b) 0 *(4,2) vs Sample Number, (0 *(4,2) = 0.96328) 
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(c) 0 *(4,3) vs Sample Number, (0 *(4,3) = -2.1772) 
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(d) 0 *(4,4) vs Sample Number, 



(0 *(4,4) = 2.3215) 
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Fig . 6 . Eigenvalues for 0 * , $ * with T * 1 , 2 , 3 = 0 , and 0 * 
using the Complete T * Resulting in Residues for the 
10 sec Dominant Period Plant Sampled and Excited at 
1 sec Intervals. No Measurement Noise. 
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SAMPLE RATE 



The choice of sample rate for discrete stochastic simulation does 

not necessarily follow the supposition that the more rapidly a system 

is sampled the "better". In some cases the existence of an optimal 

sample rate has been shown [ 6 ] . To investigate the accuracy of 

identification by batch processing, a sample set of 500 was generated 

for sample intervals from 0.02 to 4.0 seconds in increments of 0.02 

seconds. The model was excited by a different excitation set for each 

sample rate but with the same statistical properties, normal (0,1) 

th A 

"white" noise. The N row elements of 0 * for each sample set were 
subtracted from the true 0* values and plotted versus sample inter- 
val, fig. 7a,b,c,d. For sample intervals less than 0.06 seconds 
R(0) went singular. Each element appears to have its own best sam- 
ple rate ranging from 1 to 1.5 seconds for a dominant 10 second plant. 
Hence a rule of thumb: Sample at 10 times the dominant frequency! 

To eliminate the effect of varying the excitation set, the above 
procedure was repeated using the same excitation set for each sample 
period, 0.05 to 10 seconds in increments of 0.05 seconds. The re- 

A 

suiting error curves were the same as illustrated by the 0 * (4 , 2 ) curve, 
fig. 8a. The apparent stabilization of the error as the sample period 
approaches the natural period of the system is misleading. The ab- 
solute values of 0 * become small due to the large damping coefficient 
but percent errors actually increase. 

It is interesting to note that as the dominant period is changed 
the shape of the error curve is the same but shifted in the direction 
of change of the period. The ratios of sample period to natural period 
remain constant at the zero error intercept points . The error curves 

A 

for the $*( 4,2) element for natural periods of 8 and 6 seconds are 
shown in fig. 8b, c. All of these curves were generated using the full 
r * vector and hence contain numerator dynamics and their inherent 
residues. Thus it is pertinent that there is some optimal sample rate 
that minimizes residue errors. 
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Errors are not due solely to numerator dynamics. The original 
contention that the choice of sample period is significant may be em- 
phasized by repeating the above runs for the ten second plant with 
suppressed T *. The resulting identification errors are plotted versus 
sample period for each element, fig. 9a,b,c,d. Errors may be sig- 
nificantly different if the excitation, W^, drives the system at a 
different rate than the sampling (observation) frequency. 
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(b) Error = 0 *(4 , 2) - 0 * (4 , 2) vs Sample Period. 
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(c) Error = </> * (4 , 3) - d *(4 , 3) vs Sample Period. 
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(d) Error = 0 *(4,4) - 0 *(4,4) vs Sample Period. 
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Fig. 8. Error = <t> *(4,2) - £ *(4,2) vs Sample Period (0.05 to 10.0 sec) 
using Batch Processing of 500 Samples from the 10, 8, and 6 
sec Dominant Period Plants. One Sequence of 500 Discrete 
Gaussian Values, ct = 20, was used as Excitation for all Sample 
Periods and all 3 Plants. T * was not Suppressed, (a) Error of 
£ *(4,2) vs Sample Period, 10 sec Plant ^(b) Error of £*(4,2) vs 
Sample Period, 8 sec Plant, (c) Error of <b *(4,2) vs Sample 
Period, 6 sec Plant. 

(a) Error of £*(4,2) vs Sample Period (sec), 10 sec Plant. 
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(b) Error of % *(4,2) vs Sample Period (sec), 8 sec Plant. 
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§ 




(c) Error of % *(4,2) vs Sample Period (sec), 6 sec Plant. 
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Fig. 9. Error = 0 * - $ * vs Sample Period (0.05 to 10.0 sec) using 
Batch Processing of 500 Samples from the 10 sec Dominant 
Period Plant with T * Suppressed (T * 1,2,3 = 0). Discrete 
Gaussian Excitation, a = 20, was applied at the Sample Rate. 

The Same Excitation Set was used for each Different Sample 
Period, (a) Error of $ *(4,1) vs Sample Period, (b) Error of 
0 *(4,2) vs Sample Period, (c) Error of $ *(4,3) vs Sample Period, 
(d) Error of 0 *(4,4) vs Sample Period. 

(a) Error of $ *(4,1) vs Sample Period , T * 1 1 2 , 3 = 0 
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(b) Error of 0 *(4,2) vs Sample Period, r *1,2,3 
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(d) Error of 0 *(4,4) vs Sample Period , T * \ ,2 , 3 = 0 • 
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APPLICATION TO SHIP MOTION ESTIMATION 

As a result of the foregoing work and associated simulation, the 
following proposal is offered as a practical procedure for ship motion 
estimation. 

Initial identification of the order of the pitch and roll plants and 
the corresponding 0 *'s is to be evaluated for the free dynamic system 
with the ship in calm water. Initial angle may be applied with weights 
or other means. The data thus obtained may be used as a reference and 
to initialize the stochastic routines. The yaw plant may not be iden- 
tified as a free dynamic system since there is no suitable way to apply 
a controlled initial condition which will project on all the eigenvectors 
in the absence of a righting moment. 

Stochastic identification is to be applied whenever the ship is 
underway at sea. Either batch processing or the recursion equations 
may be used. Batch processing should be used at least periodically 
to verify the system order from R(0) ^ and to identify Q if required. 

The recursion equations are more suitable for time sharing and do 
not require storage to accumulate large sample sets of data. Since 
periodic identification requirements are anticipated to follow system 
dynamics , recursive identification may be initiated using the previous 

A 

0* as the best estimate and flagging in P(0) . 

The actual estimation of ship attitude follows the identification 
of the 0 * 1 s (and Q's). Equation (9) may be applied directly to predict. 
Discrete linear (Kalman) filtering is available. 

The development throughout this paper has assumed that the only 
reliable or convenient measurements available for identification were 
discrete pitch, roll, and yaw angles. If the number of observables 
is increased, the accuracy of the identifications and the estimations 
may be expected to improve [ 7 ] . 
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CONCLUSIONS 

1 . Identification of ship motion dynamics is feasible by stochastic 
methods . 

2 . The accuracy of the identification is dependent upon the numerator 
dynamics of the true system z-transfer function. The contention 
of Lee [ 4 ] that decorrelation would improve accuracy was not 
true for this work or that of Blackner [ 7 ] . 

3. The accuracy is dependent upon the sample period. An optimal 
sample rate exists. 

4 . Identification of an unexcited plant (the free dynamic system) 
having a single initial condition is independent of F* and hence 
is not affected by numerator dynamics. This should be the source 
of the most reliable identification as a reference. 

5. Further work should be done using actual ship motion data to 
complete the identification evaluation and to investigate the 
performance of discrete linear (Kalman) filtering as the final 
stage of ship motion estimation. 
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APPENDIX 



CDC FORTRAN 63 Programs and Subroutines for Digital Sim- 
ulation and Identification. 



PHITEST: 


Main program generating the simulation model 
and controlling the identification. 


PHIDENT: 


Recursive identification routine for the scalar 
observable case. 


PHICALL: 


Batch processing for general observable case. 
Identifies 0 *, N, and Q. (GAUSS3 is identical 
to RECIP) 


PHIDEL: 


Computes standard form 0 and T matrices for 
sample data systems. 


PHIDEL1 : 


Computes 0 * and T * matrices for sample 
data systems. 


PSD: 


Computes Power Spectral Density from auto- 
correlation functions . 


RECIP: 


Matrix inversion routine. 



ADD, SUB, TRANS, PROD, PRINTER, and MATREAD: 

Convenient subroutines for the indicated 



PTPLOT: 


matrix operations. 

Graphs and prints out 2 dimensional data on 
line printer. 
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THE AiOVE RELATIONSHIPS FORM THE STATE DIFFERENCE EONS 
X ( A + 1 ) - PHI*X(.<) + DEL*W(K) OR 
Y ( K.+ 1 ) = PH 1 1*Y ( K > + 0EL1*W(K) 
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PHI 1 ( I » J) = 
PHI ( I » J )=0 
H< I ,J)=0. 

A ( I » J ) =0. 
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X3 ( I + 1 ) =PH I (3»1)*X1(I)+PHI(3>?)*X2(I)+PHI(3»3)*X3(I) + 
1PHI (3>4)*X4( I )+DEL(3»l )*Wi< 

X4( 1+1 ) =PHI ( 4 , 1 ) *X 1 (.1 ) +PH I ( 4 » 2 ) *X2 ( I ) +PHI ( 4, 3 ) *X3 ( I ) + 
lPHI ( 4 > 4 ) *X4 ( I ) +DEL ( 4 » 1 )*WK 
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A 1 » 2 » 3 TEMPORARY STORAGE 
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ELEMENT OF THE PHI MATRIX 



COMPUTATION OF DEL IS COMPLETED WHEN THE TEST 
ON PHI IS SATISFIED EXCEPT FOR MULTIPLICATION BY 
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CALL PRINTER(B»N*M»ND»MD»8H B 
CALL PRINTERCPHI »N ,N » 1 0 * 10 *8H PHI 
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151 C(I*J) = C(I.J) + A( I »K)*8(K*J) 
END 



m 

co 



3 

CL 

*~ 

3 

O 

CL 

UJ 



o 

cm 



CD 

h- 

O 

u 

CL 



O 

O' 



CL 

LiJ 



CL 

CL 

LU CL 
X ID 
t- CL 
< 
Z CL 



UJ 

CD 

CL 

O 

u_ 



UJ 

X 



< 




o 




O 








— K 


























h- 






UJ 


UJ 








in 


























o 




>- 


X 










♦—4 


























-J 




-J 


1- • 


» 


































CL 




-J 


o 


UJ 








u 


























►-4 




< 


iD O 


-J 








X 






























u 


co m 


CD 


































X — 


*- 


M 


o 


< 








CL 






















rH 




• o 


O 


♦“ 


CL O 










UJ 






















vD 




>- o 


-J 


CL 


U h- 


CL 








CL 






















• 




♦ in 


0. 


UJ 


< 


< 








O 






















o 




X — 


h- 


> 


(O 


> 








CL 






















vO 




— V 


CL 




O 










a 














sO 




O' 




* 




h- * 




o 


UJ 


O 








X 


















#> 




rH 




O — 


CL 


UJ 


H- c0 


z 






















h- 




O' 




vO 




—J H 


O 


1— 


1- h- 


I— 1 








X 














•> 




♦ 




— K 




CL O 


Ll 


z 


o z 


X 






CM 


rH 














h- 




00 




X 


z 


h- m 




*— i 


-J — • 


UJ 






* 


*H 




• 










— . * 




«*■% 




< 


♦— * 


CL ~ 


<-0 


CL 


CL O 


Cl 




H 


CO 






X o X 










X 




z 




X 


X 


X 




CL 


CL 


z 






• 


X X 




O rH CM rH 










< 




►H 




>- 


V- 


UJ 


z 






l—l 


X 


CM 


CO 


* rH 




O' II -H II 


+ 






X 


X 




X 






* 


z z 


UJ 


» 


* Ll 






4 K 


— * 


o * 




«* 


X 








> 




>- 




1 


vO 


~ o 


X 


< 


UJ o 


z 


X 


CVJ 


o 


o \ 




H *-H rH »H 


rH 


rH 


rH 


rH 


1 


H 


1 


►m LU 




O UJ 


1— 


3 


U) 


h- 


< 


H 




o 


rH 




« ~ M — 


M 






M 


A 




— * 


— 3 


z 


rH 3 


3 if) 


ld 


U) 


< CL 




— «. 


CM 


in 


h- 


z 


w < w JO 


U 


>- 


>* 






>- 


*—* 


>■ Z 


•H 


Z Z 


O Z 


CL 


NH 


Z Ul 


UJ 


rH 


1 


1 


»“ < 


CL 


I- H- 


h- 


l( 


H 




'W 


H 




H *h 


X 


HOCm 


CL UJ 


< 


u 


-- CD 


CD 




X 


X 


z x 


3 ^ O m O 


O 


X 


Z O' 


>- 


X >- Z >- 


Z 3 K 


CD X 




ID 


O X 




X 




W 


•-< CL h- 


-J -J 


-J 


< 


•— • 






< 




H Z 


^■r 


H(-Z 


3 ~ 




GO 


CL 3 


X 


w 


Ll 


U. KOUJ OO.OQ. 


a. 


X 


IOI1.ZIL 


X o u. 


CL UJ O 


O) Q 




< 


o z 




X 




•— • 


CL Ll 


CL 


O O 


M 


V >- n 


•— « 


>- 


M 


>- u 


H 


CL CL U 










h- 






































II 


II H 


UJ 






rH 


CM O 




(o ^ m 












vO 




00 Ch 




O rH 






X 


V X 


-J 








O 
























vO v0 


















rH 


























U 


U 


u 


c 

c 

c 

c 

c 



































73 



PRINT 107 



ir 

c< 



oo 



X 



X 

r- 

» 

v 

x 



X 

X 



X 

lT\ 



in 

o 

H 

UJ 



(M 

X 

> 

< 



>- 

< 



< 



in 

H 

H 

•» 

4* 



0 

CM 

\ 

V 

< 

1 



X 

X 

M 



GO 

h- 

O 

-J 

a 

•— » 

u 

►- 

o 

-J 

CL 



CL 

>- 



X 



K 

CM 

m 

o 

H 

UJ 

* 

X 

CM 

m 

o 

H 

UJ 

* 

X 



CM 

CM 



* 












< 


« — . 














o 


















rH 










r— 1 






X 


o 






— - 














• 














»H 




* 








rH 


M 






o 


• 




— * 


•— « 








< 






rH 














It 




• 








+ 


►— » 






r- 


O' 




- — ■ 
















+ 


















X 
















•» 


00 




•—» 


>- 


rH 






rH 








o 












* 




rH 








<w 


«— 






in 


X 




nr 


1 








>- 








4- 
















* 








X 


»— t 




X 


• 


<r- 




>- 


Z 








< 


V 




< 


•* 












“) 




rH 








w 




— - » 


< 


o 


z 




1 


•— i 


H 








< 






o 






HD 


ao 




w> 




< 






rH 


u. 


CD 


rH 


X 


rH 


*— « 




Z 


X - 








z 






z 


•4- 






CM 


CM 




CO 




o 






CM 


lO 


►- 


< 


> 


UJ 


X 




»— « 


>- rH 


» 




O' m 


•— « 


z 










o 




* 




t— 




O' 








CO 


O co 


* 


* 


>- 




X 


— + 


rH 




CM rH 


X 


M 




X 


HD 




m 


H> CO 


O' 




o 








CM 


<M < 


-J 


H 


z 


• 


1 




>- 




H 




» •* 


V 


X 




V 


rH 




* 


cm m 


rH 




-1 




• 




CM 


CM 


1 


CL 


H 


»— t 


X 


X 






U. o 


— * 




CM cn 




>* 




w 






r- 


* * 


» 




a. 


r»» 


X 




«k 




rn 


»— t 


» 


X 


H 


< 




%r 


U) _l O 




H rH 


u. 


nr 




u_ 


O 




*-H 


rH O' O' 




»-H 


rn 


rH 




CM 


CM 


rn 




X 


V 


* 


X 




U- 


CD CL O' 




* • 


U) u. 




cO 


O' 




* 


CM rH 


rH 


CM 




w 


•> 




CM 


CM 


rn 




r— 1 


» 


X 


V 


•» 


co 


< > 


1 




CM O' 


CO </> 




CD 


1 




CM 


«•— » 




+ 




> 


X 




* 




nr 


rH 


•» 


CM 


H 




f-H CD 




I- O' 


rH CM 


< CD 




< 


VO O' 


in 


— O' 




h- CO 




rH 




O 




X 


o 


\ o 


nr 


Hr 


H 


< 


u. — 


O a> 




■W 


< 


HD 


«r> 


M 


00 




rH rH rH O O 


r»* 


w 


oo 


CM 


rH 




H 


*— 


H 




u. 


►H 


H 


a it 


-1 


n 


Z X U- 


H 


*H 


II 


X 


If 


r-* 


♦ — 


1 




*H 


Mhtn 




+ 


U- 




►- 




< 


U) 




►- 


O h Cl h *-< < 


Q 


l/) 




If) 


< lO 


1-4 


M rH 


•-4 


a. 






< 




X 


n 


UD 


h- 


< 


h- 


X 


CO 


00 O 


X O >- O X X o 


•— i 


O 


M 


X 


t—t 


w 


— 1 




>* h- 


X X o 


1 


Hr 


CD 


z 


X z 


CL 


< 


m 


-J 


— -J 


»“H 


-j 


V >- X 


X 


1— 


X 


*-H 


X 


X 


X M 


X 




z 




OC h- 


*-H 


X < 


»-H 


QC 


»— < 


O 


II 




CL 


— CL 


w 


CL 


. -w 


w 


< 




< 


w 


< 


w 


w 




M 


*-H 


— 


o 




w. 


-w< 


W 


CL O 


a: 


ll_ 


V 


O >- 


U. V lL > lL U. 


U_ 


X o 


X 


Ul X 


u. 


U. U. Ul 


X QC O 


It O li- 


U- 


Ll 


CL 


u_ 


CL 




< 


O 






























x a. 


O' 

- A 




e> 


►H 


*-H 


rn 




H 




CM 








o 




H 


rH cm cn 


<*- 




in 




o 


H) O CM CO 






C^l 




O' 


o 


rH 




O 




O 








rH 






fH rH 


rH 


rH 




rH 






h in in 


in 


rH 






o 




rH 


CM 


CM 



74 



22 IF ( I YPLOT- IX AX IS )23*24*25 

23 KK= I YPLOT + 2 
LL= I XAX IS+2 

PRINT 103 ♦( IPLOTA(J) . J = 1 . I YPLOT ) * I PLOTC » ( I PLOTA ( J ) , J*KK 
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